Modelling transmission of Mycobacterium avium subspecies paratuberculosis between Irish dairy cattle herds

Bovine paratuberculosis is an endemic disease caused by Mycobacterium avium subspecies paratuberculosis (Map). Map is mainly transmitted between herds through movement of infected but undetected animals. Our objective was to investigate the effect of observed herd characteristics on Map spread on a national scale in Ireland. Herd characteristics included herd size, number of breeding bulls introduced, number of animals purchased and sold, and number of herds the focal herd purchases from and sells to. We used these characteristics to classify herds in accordance with their probability of becoming infected and of spreading infection to other herds. A stochastic individual-based model was used to represent herd demography and Map infection dynamics of each dairy cattle herd in Ireland. Data on herd size and composition, as well as birth, death, and culling events were used to characterize herd demography. Herds were connected with each other through observed animal trade movements. Data consisted of 13 353 herds, with 4 494 768 dairy female animals, and 72 991 breeding bulls. We showed that the probability of an infected animal being introduced into the herd increases both with an increasing number of animals that enter a herd via trade and number of herds from which animals are sourced. Herds that both buy and sell a lot of animals pose the highest infection risk to other herds and could therefore play an important role in Map spread between herds. Supplementary Information The online version contains supplementary material available at 10.1186/s13567-022-01066-5.


Introduction
Bovine paratuberculosis, or Johne's disease, is a disease caused by Mycobacterium avium subspecies paratuberculosis (Map). Johne's disease is endemic in the dairy sector worldwide. It has a large economic impact due to milk losses, early culling and increased mortality [1]. Susceptibility reduces with age, and animals are usually infected in the 1 st year of life [2]. Clinical signs, including weight loss, decreased milk production, and diarrhoea, do not usually appear before first calving and are sometimes never observed [3]. However, shedding starts before animals show clinical signs [4]. Tests to identify infected animals include milk or serum ELISA, PCR, and faecal culture. However, diagnosing infected animals is difficult because of the low sensitivity of these tests, ranging from 15-71% depending on the stage of infection [5,6]. Thus, infectious animals can remain undetected in a herd. Map is mainly transmitted between herds through movement of these infected but undetected animals [7]. Understanding the effect of animal movements on the spread of Map between herds and its subsequent persistence within infected herds can be very informative for control purposes. Once Map is introduced into a herd, model Open Access *Correspondence: floor.biemans@oniris-nantes.fr Biemans et al. Veterinary Research (2022) 53:45 results suggest that a herd can remain infected for over 10 years [8,9]. Indeed, the probability of persistence 15 years after introduction has been estimated to be up to 42.7% following the introduction of a single infected animal into a typical spring calving Irish dairy herd [10].
In this study we aim to investigate Map spread between Irish dairy herds. In Ireland, the dairy industry is pasture-based to optimize the use of grass as the primary feed source for lactating cattle [11]. For more than 90% of dairy herds, cows are housed during the winter, and most calves are born in spring [11,12]. Herds have the opportunity to join the voluntary Irish Johne's Control Programme (IJCP), which has four objectives: (1) Enhance the ability of farmers to keep their herds free of Map, (2) Reduce the level of infection in infected herds, (3) Provide additional assurance to the marketplace, and (4) Improve calf health and farm biosecurity [13,14]. At the end of 2020, 11% of the dairy herds in Ireland, representing 18% of the dairy cows, were registered in the IJCP [13]. Herd prevalence was estimated to be 20.6% in 2005 [15] and 28% in 2013-2014 [16].
Although trade movements are well described in Ireland [17], the contribution and intensity of trade movements between herds on Map spread at the national scale in Ireland are still barely understood. Epidemiological modelling is relevant to answering this question, as it enables us to integrate both the within-herd population and infection dynamics, and to link herds through trade movements [18]. So far, only a few such models have been developed to predict Map spread between cattle herds at a large scale, including a model adapted to the dairy farming system of western France [19,20], a model for Northern Italy [21], and a model for Slovenia [22,23]. The French model links stochastic compartmental within-herd models to each other using cattle trade data. Recently, this model has been transformed into an individual-based model to better account for animal characteristics and to facilitate assessment of control measures [24,25]. The flexibility of this model makes it a good candidate to be adapted to represent Map spread in the predominantly seasonal Irish dairy farming system [26].
We used an adapted version of the French transmission model to simulate Map spread. Our objective was to investigate the effect of herd characteristics on Map spread in dairy herds in Ireland. Herd characteristics that could be calculated from the movement data included herd size, number of male animals introduced for breeding, number of animals purchased and sold, and number of herds that the focal herd purchases from and sells to. We used these characteristics to classify herds in accordance with their probability of becoming infected and of spreading infection to other herds.

Materials and methods
A stochastic individual-based model was used to represent herd demography and Map infection dynamics of each dairy cattle herd in Ireland [10,24]. Data on herd size and composition, as well as birth, death, and culling events were used to characterize herd demography. Herds were connected with each other through observed animal trade movements [17]. A comprehensive dataset on cattle trade movements at a national scale in Ireland was used, comprising information for each moved animal (e.g., age at move, sex, breed, date of move, reason for the move, and source and destination herd). Details on the datasets used, and on the processes occurring at withinand between-herd scale are presented hereafter.

Data
Data were extracted from the Animal Identification and Movement system (AIM) of the Irish Government's Department of Agriculture, Food and the Marine. The AIM database comprises records on bovine births, movements and movement types [17,27]. Movement data and demographic data (e.g., animal ID, sex, breed, date of birth) were extracted for the 10-year period 1 st January 2009-31 st December 2018.
Herds included in our dataset, also referred to as the studied metapopulation, were dairy herds classified as such using a similar methodology to Brock et al. [28]. In general, European dairy herds are organized into groups of animals that are a similar age. Each group has its own housing and management that influence how, and with which other animals, they can interact [29]. The model was developed to simulate this contact structure for typical European dairy herds [9,24,30]. Therefore, only dairy herds were included.
Only herds with more than half of their animals of dairy breed (Additional file 1) were selected. Furthermore, more than 30% of the animals had to be female. With these two criteria, we were able to separate mixed herds (with female and male animals of dairy and beef breed) from store dairy male herds (high proportion of males that are of dairy breed [28]) and from store beef mixed herds (high proportion of beef animals that could be either males or females [28]). Next, herds were selected based on completeness of movement data (data were considered complete when data were available for a herd for each year of the study period). To simplify the simulations and analyses, we chose to exclude herds that closed or started during the study period (2127 herds). A herd was considered to be closed when, over the entire study period, animals only entered the herd via birth and exited the herd via death or slaughter. Because these herds do not trade with other herds, they do not contribute to Map spread and were excluded from the dataset.
Only female animals were modelled, with the exception of the purchase of some male animals for breeding (see below). Hence, herds that did not exchange female animals with other herds were also excluded. Finally, herds were selected on the basis of size: herds with 15 animals or less, or five adults or less, were considered to be non-commercial dairy herds and were excluded. A total of 13 353 herds, with 4 494 768 dairy female animals, were included in the metapopulation. The retained herds belonged to one of six possible herd types (defined in Table 1), as described in Brock et al. [28]. Figure 1 presents the distribution of herds included in the metapopulation per herd type. Some farms within the metapopulation conducted both dairy and fattening activities, however, for the purpose of modelling, it was assumed that the fattening activities took place on a different part of the farm with negligible contact between the dairy and fattening herds, and thus the fattening activities were not represented in the model.
In Ireland, the use of natural service bulls following rounds of artificial insemination is common practice [31]. To account for the risk of introducing Map into a herd via the purchase of bulls for breeding purposes, some males were also included in the dataset. Including bulls was only considered for herds that purchased dairy or beef bulls (83% of the herds). To differentiate bulls introduced for fattening purposes from bulls introduced for breeding purposes, only bulls that met the following criteria were included: bulls had to be pure bred dairy or beef, they had to enter the herd via trade, and beef bulls had to come from a non-dairy breeder herd. Bulls had to be at least 1 year of age upon entering a herd, and they had to remain in a herd for at least 5 months (one breeding season). The maximum length of stay in a herd was 4 years for beef bulls and 2 years for dairy bulls. It was assumed that dairy bulls were only used for breeding for a maximum of 2 years because after that these animals had the potential to mate with their daughters (who first calve at approximately 2 years of age). The number of bulls introduced per year was calculated based on the assumption that about half the cows would be serviced by a bull, with 20 cows per bull. The other half of the cows will be artificially inseminated. The number of cows present in a herd was determined on the first of May of a given year. Based on this calculation, if more bulls were needed than present in the trade data, dairy bulls that were born into a herd and kept until breeding age were selected. A total of 72 991 bulls were selected for 11 121 herds of the metapopulation and added to the dataset.

Within-herd Map transmission model
The within-herd model is adapted from Camanes et al. [24] and Biemans et al. [10]. All parameters used and model equations related to disease transmission are described in detail in the supplementary material of Table 1 List of herd types A full description of these herd types can be found in [28].  [10]. Please refer to that study for a detailed overview. Briefly, it is a stochastic individual-based model with a discrete time step of 1 week. A time step of 1 week was chosen because it makes it possible to accurately represent the susceptible decay and the seasonal herd management of Irish dairy herds while reducing computation time compared to a time step of, e.g., 1 day. It accounts for herd structure and infection dynamics (Additional file 2). Animals belong to one of six age groups: newborn calves, unweaned calves, weaned calves, young heifers, bred heifers, and cows. Depending on the age group that animals are in, they move to the next age group either because they reach a defined age or at a defined time in the year (Additional file 2). Animals also belong to one of six health states: susceptible (S), resistant (R), transiently infectious (I T ), latently infected (I L ), moderately infectious (I M ), and highly infectious and with the possibility to also be clinically affected (I H ). Animals are assumed to be most susceptible at birth, with susceptibility decreasing exponentially with age [2]. Infectious animals (I L , I M , I H ) shed Map in their faeces and, after calving, in their colostrum and milk. The quantity shed depends on the health state and is heterogeneous between animals of the same state [4,32]. Animals can get infected via the following transmission routes: in utero including transmission during parturition, via ingestion of contaminated colostrum or milk (directly or indirectly with faeces), via contact with the local environment or the general indoor environment contaminated with faeces. Transmission via the local environment is defined as the risk posed by other animals held in the same place but not necessarily at the same time. The local environment can be indoor or on pasture. Transmission via the general indoor environment is defined as the risk posed by other animals held indoors but not necessarily at the same place or time. All infectious animals that reside indoors contribute to the contamination of the general indoor environment. Within the indoor environment, 40% of the Map load present was removed per week, representing the effect of removing manure [9]. Similarly, on pasture, a 7.1% reduction in Map load occurred weekly [33]. Additionally, when there are no animals present in one of the indoor environments, they are cleaned more thoroughly, removing an extra 16.7% of the Map load [10,30]. Details of all model parameters are in Additional file 1.

Herd type Abbreviation Description
To represent Irish dairy farm management, which is highly seasonal [26], animals are divided into two clusters depending on the season in which they are born in order to prevent mixing of yearlings. The spring cluster consists of animals born in the first half of the year, while the autumn cluster consists of animals born in the second half of the year. The cluster, animal age, and the time of the year determine the age group and environment of an animal, e.g., spring-born unweaned calves are kept indoors up to week 14 of the year and on pasture from week 14 onwards, whereas autumn-born unweaned calves are always kept indoors. For both clusters, details of the transitions between age classes and environments in which animals are kept are in Additional file 2.
For each herd in the metapopulation, herd parameters for size and exit rates were calculated from the data. Herd size was calibrated on the 1 st of January 2009. For each age class, seasonal exit rates (i.e., exit rates over a 3-month period; January-March, April-June, July-September, October-December) were defined because it resulted in the best agreement between observed and modelled herd size compared to the use of yearly or monthly exit rates (Additional file 3). A newborn calf was added to the herd whenever a birth was observed in the data.

Between-herd Map transmission model
The between-herd model connects all the within-herd dynamics through animal movements. Movements from and to a herd in the metapopulation are modelled as observed in the data. When an animal was moved via a market, the movement was represented as if it occurred directly from the source herd to the destination herd. The movement dataset comprises 2 304 149 animal movements of which 17% were between herds in the metapopulation, 20% were from a herd outside of the metapopulation to a herd in the metapopulation, and 63% were from a herd in the metapopulation to a herd outside of the metapopulation.
The movement data define the date that a movement occurs, the age of the animal moved, and the source and destination herd. The animal to be moved is randomly selected from the relevant age group in the source herd, and thus can be of any health state. When there is, due to chance, no animal present in the correct age group, an animal of the closest age group is selected. If an animal is coming from a herd outside of the metapopulation, its health state is drawn from a distribution that corresponds to the proportion of animals in each health state in the same age group within the entire metapopulation. Thus, it is assumed that the average risk of introducing an infected animal is the same within and from outside of the metapopulation.

Herd characteristics
Seven herd characteristics are calculated from the movement data for every year: in-degree, out-degree, instrength, out-strength, polarity, herd size, and number of male animals introduced. The in-degree measures the number of herds that the herd of interest receives animals from. The out-degree measures the number of herds that the herd of interest sends animals to. The instrength measures the number of animals purchased by the herd of interest (incoming movements). The outstrength measures the number of animals sold by the herd of interest (outgoing movements). Polarity is defined as Polarity = #incoming animals−#outgoing animals #incoming animals+#outgoing animals . Polarity takes values between −1 and 1 and represents the trading behaviour of a herd, where herds with a polarity between −1 and −0.25 can be considered as predominantly selling herds, between −0.25 and 0.25 as neither predominantly buying or selling herds, and between 0.25 and 1 as predominantly buying herds [34]. The herd size was calculated as the total number of animals present in a herd on the first of January. These herd characteristics are calculated for every year and then averaged over the 10-year period. The annual strength and degree can vary over the years. The median values for the difference between the lowest and the highest value were 2.0 for in-degree and out-degree, 9.0 for in-strength, and 18.0 for out-strength (Additional file 4).

Simulation settings
Which herds were chosen to be infected at the start of the simulations had an effect on model predictions (Additional file 5). Therefore, we first assessed which herds were the most likely to be infected 10 years after Map introduction in the metapopulation when they were not the ones initially infected. On this basis, for 1000 replicates, 25% of the herds were randomly chosen to be initially infected irrespective of their herd characteristics. The within-herd prevalence is low in the majority of herds based on earlier analyses of field data [15,35]. Therefore, the within-herd prevalence in the herds that are initially infected was drawn from a Gaussian distribution N(−0.42,0.12), keeping only values below 0.7 because higher values are rarely observed in the field (Additional file 6). The probability of a herd being infected after 10 years of simulating Map transmission was calculated for each herd given that it was not initially infected. The 30% most likely herds to be infected were chosen as candidate herds to be initially seeded with Map in the subsequent simulations.
Second, for 300 replicates, amongst the candidate herds that were most likely to be infected, 25% of the total number of herds were chosen to be initially infected. Their initial within-herd prevalence was drawn from the aforementioned distribution. Assuming a herd prevalence of 25% in 2009 was in agreement with field observations, the prevalence being estimated at 20.6% in 2005 [15], while it has been estimated at 28% in 2013-2014 [16]. Map transmission was simulated for 10 years, matching the temporal extent of the movement data.

Model outputs
Firstly, we analysed over time the proportion of herds infected, hereafter called "herd prevalence", and the proportion of infected animals among > 2-year-old animals present in each infected herd the first of January each year, hereafter called "within-herd prevalence". Secondly, we investigated the correlation between observed herd characteristics to provide context when interpreting the results of the simulations. Thirdly, we calculated for each herd the probability of becoming infected and the probability of escaping infection. We related these probabilities to the in-and out-degrees, in-and out-strengths, herd size, and number of male animals introduced for breeding. The probability of a herd becoming infected was calculated with as denominator the number of replicates in which a herd was not initially infected and as numerator the subset of these replicates in which at least one infected animal entered the herd during the 10-year simulation period. The probability of a herd escaping infection was calculated using the same denominator (the number of replicates in which a herd was not initially infected) and as numerator the subset of these replicates in which no infected animal entered the herd within the 10-year simulation period. A generalized linear model with a logit link was used to connect (the logarithm of ) each herd characteristic to the probability of a herd being infected, model fit was assessed by Akaike information criterion (AIC). Fourthly, we investigated the infection risk posed by each selling herd type. We determined the number of buying herds infected following outward movements of infected animals from each selling herd, and summarised as the average number of herds infected by each selling herd type. A buying herd was considered to be infected by the selling herd if at least one infected animal was introduced and up to that moment no infected animals had been present in the herd. Lastly, we determined the number of unique infection sources per herd. A herd is an infection source if it introduces an infected animal in the herd of interest in at least one replicate.

Prevalence
For all replicates, herd prevalence increased over time (Figure 2A). After 10 years of simulation, herd prevalence was 49.9% on average. The average within-herd prevalence among > 2-year-old animals for all herds (including uninfected herds) was 4.1% at the start of the simulations in January 2009. It increased to 5.6% at the end of the simulations in December of 2018. The average withinherd prevalence among > 2-year-old animals within infected herds ( Figure 2B) was 17.8% in 2009, decreasing to 8.5% in 2012 then increasing to 12.7% at the end of 2018. This initial decrease was expected as only few infected animals are present in newly-infected herds, which by calculation induces a decrease in the average within-herd prevalence. Table 2 presents the mean value of each herd characteristic per herd type. Herds classified as SdR (herds that rear dairy females) have the highest values for in-strength and out-strength but the number of SdR herds is small (n = 2). On average, herds classified as DnR-C (dairy herds which contract out calf rearing to other herds) have a higher out-degree and out-strength

Table 2 Mean value of a herd characteristic per herd type
The 25% and 75% quantile values are in parentheses. Herd types are typical dairy (D), dairy no rearing-contract (DnR-C), dairy no rearing-no contract (DnR-nC), dairy also rearing male calves (DRm), mixed (M), and herds that rear dairy females (SdR) (defined in Table 1).

Herd type
In compared to the other herd types. Furthermore, herds classified as DnR-C or DnR-nC (non-rearing dairy herds which do not use contract rearing) have a higher in-strength compared to the other herd types. Figure 3 presents the correlation between the six herd characteristics for all 13 353 herds. In-degree and instrength have a strong positive correlation (0.80), as have in-strength and out-strength (0.72). For the other characteristics, the correlation is moderate (0.36-0.67) to weak (< 0.36). Figure 4 presents the probability of a herd becoming infected versus six herd characteristics. For all of these herd characteristics, the probability to becoming infected increased with higher values for the characteristics. For all models, the explanatory variables were significant, but the model with the logarithm of in-strength as the explanatory variable had the best fit (lowest AIC). For 957 herds, the probability to becoming infected was zero. However, this did not necessarily mean that they did not purchase any animals, as 177 of these herds had an average in-degree (mean = 0.15) and in-strength (mean = 0.36) higher than 0. For 2390 herds, the probability to becoming infected was 100%. However, this did not necessarily mean that they purchased animals from many other herds, as 481 of these herds (20.1%) had an average in-degree lower than 1 (mean = 0.69). Herds with an average in-degree higher than 3, herds with an in-strength higher than 8 and herds that purchased more than three breeding bulls per year all had a probability exceeding 95% of becoming infected at least once during the 10 years of simulation.

Probability of becoming infected and probability of escaping infection
Of the 13 353 herds in the metapopulation, 1969 herds escaped infection in more than 90% of the replicates. Figure 5 presents the distribution of herd characteristics among these 1969 herds compared to all herds in the metapopulation, and to herds that became infected during the 10-year simulation period in every replicate. For every herd characteristic, the average value was the lowest for herds that escaped infection and the highest for herds that always became infected. This difference was most pronounced for the average in-strength, with a median value of 0.1 for herds that escaped infection and 14.0 for herds that always became infected. Of the herds that escaped infection, 57.2% never purchased a male animal, while only 4.4% of the herds that were always infected never purchased a male animal. For 99.5% of the herds that escaped infection, the average in-strength was lower than 4, while for 96.2% of the herds that always became infected the average in-strength was 4 or more. Figure 6 presents the distribution of the infection risk posed by each selling herd type. Herds classified as DnR-C infected significantly more dairy herds on average compared to all other herd types. Herds classified as DnR-nC and herds classified as typical dairy (D) infected significantly more herds than mixed herds (M) and dairy herds that also rear male calves (DRm). DRm herds infected significantly more dairy herds than mixed herds. Thus, for the mean average number of dairy herds infected per herd type: DnR-C > DnR-nC and D > DRm > M. The percentage of herds per type that did not infect any other dairy herds was 23.0% for DnR-C, 48.2% for DnR-nC herds, 54.1% for D, 65.6% for DRm, and 74.1% for M. The DnR-C herds that infected other dairy herds infected 1.7 dairy herds on average, with a maximum of 29.9. For all herd types except SdR, herds that infected other dairy herds infected between 1.1 and 1.3 herds on average. The mean of herds classified as SdR (n = 2) was not significantly different from any other herd type. Figure 7 presents the average number of susceptible herds in which the herd of interest (the selling herd) introduced an infected animal versus the average polarity of the herd of interest. Herds with an average in-strength of less than 4 are presented separately from herds with an average in-strength of 4 or more, based on the results presented in Figure 5C. Herds could have a very different in-strength while having the same polarity, e.g., a herd that sells one animal and buys one animal has a polarity of 0, while a herd that sells 50 animals and buys 50 animals also has a polarity of 0. However, the risk of becoming infected and spreading infection is not the same. The mean polarity of herds with an in-strength of less than 4 was −0.49, and for herds with an in-strength of 4 or more it was −0.09. The mean probability of becoming infected during the 10-year simulation period was 40.0% for herds with an in-strength of less than 4, and 93.2% for herds with an in-strength of 4 or more. The mean number of susceptible herds into which an infected animal was introduced was 0.34 for selling herds with an in-strength of less than 4 and 70% of these selling herds did not infect any other herd. For selling herds with an in-strength of 4 or more this number was 0.89, and 41% of these selling herds did not infect any other herds. Figure 8 presents the distribution of the number of unique infection sources per herd type. DnR-nC herds had significantly more unique infection sources compared to the other herd types. M herds had significantly more unique infection sources than D, DRm, and DnR-C. DRm herds had significantly more unique infection sources than D and DnR-C, and D herds had significantly more infection sources than DnR-C. So, for the mean number of unique infection sources: DnR-nC > M > DRm > D > DnR-C. Herds had on average 3.3 unique sources of infection, with a maximum of 20. Again the mean of herds classified as SdR was not significantly different from other herd types.

Discussion
We used a stochastic individual-based transmission model to simulate Map spread between dairy cattle herds in Ireland, with the aim of investigating the effect of herd characteristics on the risk of becoming infected and spreading infection on a national scale. We showed that the probability of an infected animal being introduced into the herd increases both with increasing in-strength (the number of animals that enter a herd via trade) and in-degree (the number of herds from which animals are sourced). Furthermore, herds that both buy and sell a lot of animals pose the highest infection risk. Non-rearing dairy herds trade more animals and with more dairy herds than other herd types and could therefore play an important role in Map spread between dairy cattle herds.
We found that herds classified as non-rearing dairy herds infected on average more herds compared to other herd types. These herds (DnR-C and DnR-nC) do not rear their own female dairy calves. Calves are moved to contract rearing herds and pregnant heifers are bought back from the same contract rearing herd (DnR-C). Or cows are bred to beef bulls and the calves are sold, replacement heifers are bought from other herds (DnR-nC). The instrength and out-strength (and to a lesser extent in-degree and out-degree) of these herd types was high compared to other herd types. DnR-C and DnR-nC herds also infected on average more herds compared to other herd types, and DnR-nC herds had the highest number of unique infection sources. Over 84% of the DnR-C and over 52% of the DnR-nC herds infected at least one other herd after they became infected themselves. Furthermore, 89% of the DnR-C herds and 40% of the DnR-nC herds had a probability greater than 97.5% of becoming infected during the 10-year simulation period. Because of these herd characteristics, it could be that these herd types play an important role in disease spread between herds. Although the trading structures are very different, pig herds in the pork supply chain with a high out-degree or long outgoing infection chain were of particular importance during an epidemic with regard to disease spread [45]. The in-going and out-going infection chains, while being complex to calculate for a national trade network of the size of the one considered in our study, could provide interesting insight in how much every herd type actually contributes to Map spread between dairy cattle herds.
We observed a strong positive correlation between indegree and in-strength (0.80), meaning that herds that bought from a small number of other herds were also likely to buy a small total number of animals and vice-versa. Furthermore, we found a high correlation between in-strength and out-strength (0.72), meaning that herds that bought a lot of animals were likely to also sell a lot of animals, and vice-versa. Tratalos et al. [17], who analysed movements of all cattle in Ireland in 2016 including movements of males and animals of beef breed, also observed a strong correlation between in-degree and in-strength (0.88). However, there was little correlation between in-and out-strength (0.057). This could be because the trading pattern of female animals for dairy herds is different from the trading pattern of all animals (including males and animals of beef breed) for all cattle herds in Ireland which involves many different herd types. For most dairy herds trade is from herd to herd. These movements are all considered when calculating the in-and out-strength. In contrast, many beef herds sell directly to abattoirs, these movements were not considered for the calculations [17].
We found that the probability of a herd becoming infected during the 10-year simulation period increased rapidly with increasing number of animals purchased per year. Herds that bought more than 8 animals per year had, on average, a probability of over 95% of becoming infected. Also for France the probability of becoming infected increased with the number of animals purchased; French farms that bought three or more animals per year had a probability of over 50% of becoming infected within the 9-year simulation period [19]. However, buying an infected animal does not always imply disease persistence in the herd. For France, the probability of persistence decreases when the number of animals traded increases, i.e., a high turnover rate increases the probability of removing infected animals [19]. Further work could focus on investigating whether this phenomenon is also observed in Ireland.
We used polarity as a measure to investigate the probability of a herd becoming infected and thereafter being a risk of spreading infection to other herds. Results are presented separately for herds with an average in-strength of less than 4 and herds with an average in-strength of 4 or more. Herds with an average in-strength of less than 4 often escape infection and herds with an average in-strength of 4 or more have a high probability to become infected (Figures 4C and 5C). There is also a difference in mean polarity for herds with an in-strength of less than 4 and herds with an in-strength of 4 or more. Herds with an in-strength of less than 4 are generally sellers, while herds with an instrength of 4 or more are generally wholesalers (herds that both buy and sell animals). For France, it was shown that wholesalers are more likely to spread infection compared to sellers or buyers [19]. In Ireland this is especially true for wholesalers with a high in-strength; because of the large number of animals going in and out of these herds, they not only have a high probability of becoming infected, but also have a high probability of spreading infection to other herds ( Figure 7).
Map infection was very persistent, as shown by herd prevalence increasing over time for all replicates. Model predictions in dairy herds in France showed that, even with low herd prevalence and within-herd prevalence, extinction was impossible without intervention strategies [19]. However, herd prevalence could be lowered using combinations of control measures, for example by combining intervention strategies with risk based trading based on herd status to prevent Map introductions in free herds [20,25]. In Ireland, herd prevalence was estimated at 28% in 2013-2014 [16]. In the French modelling study looking at risk-based trading [25], a much higher initial herd prevalence was assumed, but considering very low within-herd prevalence, most infected herds had a within-herd prevalence below the detection threshold. Another major difference between the two modelled regions lies in the number of trade movements, with many more trade movements in Ireland (2 304 149 for 13 353 herds) than in France (919 304 for 12 857 herds). It would, therefore, be interesting to investigate whether combining risk-based trading with intervention strategies would lead to similar results to those for France.
We included breeding bulls in the dataset because we wanted to account for the associated risk of Map introduction. After their introduction in the herd, bulls were assumed to behave like, and reside in the same location  as, female animals of their age group. Transmission via semen [46][47][48], whether through natural service or artificial insemination, was not modelled as such. Beef bulls that were introduced came from herds outside of the metapopulation. The probability of a herd becoming infected increased with an increasing average number of bulls purchased per year. However, the average number of bulls purchased per year was not the best predictor for the probability of a herd to become infected, average instrength being better. The number of bulls purchased per year was low (mean = 0.56) compared to the in-strength (mean = 6.37). Therefore, to reduce the risk of becoming infected, it is probably more effective to lower the number of female animals introduced from other herds than to stop using beef bulls for breeding.
The model is data-driven, meaning that herd-specific parameters such as initial herd size, exit rate, calves born over time and animal movements between herds were based on real trade data. An advantage of using real data is that there is a close resemblance to reality; actual herd demographics resemble simulated herd demographics, and herds interacting in reality are also interacting in the model. It was, therefore, possible to identify which herds with associated herd characteristics were most exposed to infection, had a higher chance of escaping it, or contributed most to Map spread. A disadvantage of using real data is that the simulated period and the number of herds modelled are constrained by the data, and simulations are based on past events. Farm numbers declined over the studied period [49,50], and only herds with data over the whole period 2009-2018 were included in the model. Predicting animal movements based on trade patterns between herds would make it possible to simulate future dynamics as well. However, no method is available so far to predict animal movements at such a large scale (several thousand herds). In addition, obtaining accurate predictions might be difficult for some herds that have highly variable trade characteristics between years.
In the model, several assumptions are made. Validating the realism of these assumptions is not always easy because detailed data is not available. However, model components and the parameters used in the model are based on research, and updated when new information becomes available. For the within-herd part of the model, the transmission rates are the most uncertain. Two previous studies using a similar model, performed a sensitivity analysis on these parameters. For transmission within Irish dairy herds, only variation in the transmission rate parameter for the general indoor environment had an effect on prevalence over time, but the conclusion that the general indoor environment was the most important transmission route remained unchanged [10]. Similar results were found for transmission within French dairy herds, when transmission rate parameters were varied with 50% compared to the reference scenario [9]. For the between-herd part of the model, the only assumption that was made is that the herd and within-herd prevalence of herds outside of the metapopulation is the same as the herd and within-herd prevalence of herds within the metapopulation. Data to validate this assumption are lacking, and therefore it was simplest to assume that the average risk of introducing an infected animal from within and from outside of the metapopulation was the same.
Brock et al. [28] identified 17 different herd types in Ireland based on nine variables. Our focus was on dairy herds, and therefore, only herds belonging to one of the following six dairy herd types were included in the model. In typical dairy herds (D), the most common dairy herd type, most female dairy calves are reared to become replacement heifers and most of their male calves are sold at an early age. Therefore, there are almost no males between the age of 1 and 2 years in these herds. Mixed herds (M), the second most common dairy herd type, are characterised by having both milk and beef production activities; on average, these herd have half purebred dairy females and half cross-bred dairy and beef animals. For these mixed herds, even though there are two production types on the same farm, we assumed that the milk and beef production activities were managed in separate locations and that contact between dairy and beef animals was negligible. Therefore, we also neglected in-and out-movements related to the beef production activity. The third most common herd type is dairy herds that also rear their male calves (DRm). Similar to typical dairy herds (D), female dairy calves are reared to become replacement heifers, however, in DRm herds the males are kept, resulting in a high proportion of male animals between the age of 1 and 2 years. In the model, only data on female animals was included (with the exception of some data on breeding bulls), therefore, data on male animals in DRm herds was not included. Since the males are sold as young stock before the age of 104 weeks for fattening [51], we assumed that the probability that they are shedding Map substantially in the presence of young dairy female calves is low, and therefore, their contribution to Map spread is minimal. Slightly less prevalent were the non-rearing dairy herds (DnR-C and DnR-nC). These herds sell most of their calves, with female dairy calves being moved to external contract rearing herds (e.g., herds that rear dairy females (SdR)). The difference between DnR-C herds and DnR-nC herds is that in DnR-C herds the female calves return to their birth herd as pregnant heifers (contract rearing), whereas in DnR-nC herds there is no contract rearing and replacement animals are bought from other herds with a surplus of cows or pregnant heifers. In SdR herds, female dairy calves are reared and inseminated before returning to their birth herd (the DnR-C herds). In these herds, animals are usually of young age. Since herds with less than five adults were too small to be taken into account by the model, only two SdR herds were included in the metapopulation.
Non-dairy herds were connected to the dairy herds as well. However, 63% of the movements were from a dairy herd in the metapopulation to a herd outside of the metapopulation. Only 20% of the movements were from a herd outside of the metapopulation to a herd in the metapopulation. Brock et al. [51] investigated transport flows per herd type for Ireland. The majority of the movements to dairy herds came from store herds. In our manuscript, only two store herds (SdR) were included, mainly because in store herds often only young animals are present (i.e., young female dairy calves are introduced, reared and inseminated, before being returned to their birth herd). Most store herds were excluded from the metapopulation based on the criteria that at least five adult animals needed to be present in the herd. However, movements from these herds into the metapopulation were still simulated. The risk of introducing an infected animal through such a movement was drawn from a distribution that corresponded to the proportion of animals in each health state in the same age group within the entire metapopulation. So the risk of transmission was accounted for. Unfortunately, it was not possible to include non-dairy herd types because the model was specifically designed to represent dairy herd management.
In the model, we assumed that herd management was similar for all farms, even though the characteristics of the herd types, especially with regard to which age groups are present, are different. In the model, animals were divided into age classes, each with their own properties, and therefore all herds could be modelled in the same way. For example, for all herds we assumed that cows would go to pasture in the beginning of March, whereas calves stay indoors to the end of April. Thus, properties between age classes were different but for a given age class properties were the same across farm types.
We have shown that herds that buy many animals from a lot of different source herds have the highest probability of becoming infected, whereas herds that buy only a few animals from a few sources have the highest probability of escaping infection. Furthermore, on average, herds classified as non-rearing dairy herds are at greatest risk of infecting other herds, compared to other herd types. These findings could be used to identify herds on which the control programme could particularly focus, including those at highest risk of being infected and those at highest risk of spreading infection to other herds. For example, by combining intervention strategies with risk based trading, there is the possibility to prevent Map introduction into herds that are currently free of infection.